import geopandas as gp
import folium
from shapely import geometry as gm
import itertools
from polygon_geohasher.polygon_geohasher import polygon_to_geohashes, geohashes_to_polygon
import pandas as pd
pd.set_option('display.max_columns',None)
Public land datasets can be found at https://www.sciencebase.gov/catalog/item/62226321d34ee0c6b38b6be3 which comes from the PAD-US files, you can learn more about PAD-US dataset here, https://www.usgs.gov/programs/gap-analysis-project/science/pad-us-data-overview.
The roads dataset comes from USGS and The National Map. The dataset is here https://prd-tnm.s3.amazonaws.com/StagedProducts/Tran/Shape/TRAN_Wyoming_State_Shape.zip and you can learn more about the transportation dataset here https://www.sciencebase.gov/catalog/item/5a61c942e4b06e28e9c3bddc
The motivation to provide this analysis stems from a recent corner crossing case (https://wyofile.com/ranch-owner-corner-crossing-damages-could-exceed-7m/). HuntScore is providing this analysis and code as a means to help promote solutions to making public land accessible. In particular, an incosequential amount of public land is either entirely surrounded by private land or requires access via corner hopping (a still legal gray area). Hopefully this analysis and code can be used as a starting point to prioritize access through either public or private initiative.
Hopefully this notebook serves as an example of how to identify different access modes by using graph networks. In addition, the graph structure could be useful to weight potential access paths, where more than one exists, to provide access to the optimal amount of public land.
While the legality of corner crossing works its way through Federal court, there could many uses depending on the outcome. If found legal with a narrow prescriptive ruling where persons crossing a corner must be certain of the corner, the analysis could be used to identify the largest tracts of public land that could be accessible through corner crossing and help prioritize the placement of visible survey corners along the access paths. On other hand, if corner crossing is found to be illegal unless consent is obtained by all bordering property owners, a private or public entity could work to negotiate access with private property owners corners by merging in a private landowner dataset to this one.
wyland = gp.read_file('/Users/mhabiger/Downloads/PADUS3_0_State_WY_SHP/PADUS3_0Fee_StateWY.shp')
wyland.head(3)
| FeatClass | d_FeatClas | Category | d_Category | Own_Type | d_Own_Type | Own_Name | d_Own_Name | Loc_Own | Mang_Type | d_Mang_Typ | Mang_Name | d_Mang_Nam | Loc_Mang | Des_Tp | d_Des_Tp | Loc_Ds | Unit_Nm | Loc_Nm | State_Nm | d_State_Nm | Agg_Src | GIS_Src | Src_Date | GIS_Acres | Source_PAI | WDPA_Cd | Pub_Access | d_Pub_Acce | Access_Src | Access_Dt | GAP_Sts | d_GAP_Sts | GAPCdSrc | GAPCdDt | IUCN_Cat | d_IUCN_Cat | IUCNCtSrc | IUCNCtDt | Date_Est | Comments | SHAPE_Leng | SHAPE_Area | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Fee | Fee | Fee | Fee | FED | Federal | BLM | Bureau of Land Management | None | FED | Federal | BLM | Bureau of Land Management | BLM | PUB | National Public Lands | None | Billings Field Office | Billings Field Office | WY | Wyoming | USGS_PADUS3_0Fee_BLM_SMA_ADMU_Union | SMA_for_PADUS_20210120_1.gdb\SurfaceManagement... | 2021/06/11 | 1 | None | 0 | OA | Open Access | GAP - Default | 2021 | 3 | 3 - managed for multiple uses - subject to ext... | GAP - Default | 2021 | Other Conservation Area | Other Conservation Area | GAP - Default | 2021 | None | None | 1.688651e+04 | 2.326939e+03 | MULTIPOLYGON (((-964334.087 2508139.311, -9643... |
| 1 | Fee | Fee | Fee | Fee | FED | Federal | BLM | Bureau of Land Management | None | FED | Federal | BLM | Bureau of Land Management | BLM | PUB | National Public Lands | None | Billings Field Office | Billings Field Office | WY | Wyoming | USGS_PADUS3_0Fee_BLM_SMA_ADMU_Union | SMA_for_PADUS_20210120_1.gdb\SurfaceManagement... | 2021/06/11 | 0 | None | 0 | OA | Open Access | GAP - Default | 2021 | 3 | 3 - managed for multiple uses - subject to ext... | GAP - Default | 2021 | Other Conservation Area | Other Conservation Area | GAP - Default | 2021 | None | None | 7.044126e+03 | 8.331714e+02 | MULTIPOLYGON (((-864957.100 2495472.292, -8649... |
| 2 | Fee | Fee | Fee | Fee | FED | Federal | BLM | Bureau of Land Management | None | FED | Federal | BLM | Bureau of Land Management | BLM | PUB | National Public Lands | None | Buffalo Field Office | Buffalo Field Office | WY | Wyoming | USGS_PADUS3_0Fee_BLM_SMA_ADMU_Union | SMA_for_PADUS_20210120_1.gdb\SurfaceManagement... | 2021/06/11 | 867546 | None | 0 | OA | Open Access | GAP - Default | 2021 | 3 | 3 - managed for multiple uses - subject to ext... | GAP - Default | 2021 | Other Conservation Area | Other Conservation Area | GAP - Default | 2021 | None | None | 8.825345e+06 | 3.510835e+09 | MULTIPOLYGON (((-733353.658 2316434.437, -7333... |
wyland.groupby('d_Own_Name').GIS_Acres.sum().sort_values(ascending=False)
d_Own_Name Bureau of Land Management 17804113 Forest Service 9226163 State Land Board 3561307 National Park Service 2330195 Bureau of Reclamation 904445 State Fish and Wildlife 388974 U.S. Fish and Wildlife Service 65271 Non-Governmental Organization 31197 State Park and Recreation 18593 Other or Unknown State Land 14877 City Land 12548 County Land 7548 Private 1105 Other or Unknown Local Government 420 Regional Agency Land 389 Joint 9 Name: GIS_Acres, dtype: int64
wyland = wyland[['Own_Name','geometry']]
wyland['geometry1'] = wyland.geometry.to_crs({'init': 'epsg:4269'})
wyland['geometry'] = wyland.geometry
/Users/mhabiger/anaconda3/envs/basic_huntscore/lib/python3.8/site-packages/pyproj/crs/crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string)
wyland.geometry.values[3]
osm = folium.Map([42.994,-105.590],zoom_start=8)
folium.GeoJson(gm.mapping(wyland.geometry1.values[3])).add_to(osm)
osm
wyland.crs
<Derived Projected CRS: ESRI:102039> Name: USA_Contiguous_Albers_Equal_Area_Conic_USGS_version Axis Info [cartesian]: - [east]: Easting (metre) - [north]: Northing (metre) Area of Use: - undefined Coordinate Operation: - name: unnamed - method: Albers Equal Area Datum: North American Datum 1983 - Ellipsoid: GRS 1980 - Prime Meridian: Greenwich
wyroads1 = gp.read_file('/Users/mhabiger/Downloads/TRAN_Wyoming_State_Shape/Shape/Trans_RoadSegment.shp')
wyroads2 = gp.read_file('/Users/mhabiger/Downloads/TRAN_Wyoming_State_Shape/Shape/Trans_TrailSegment.shp')
According to technical documentation, https://www2.census.gov/geo/pdfs/maps-data/data/tiger/tgrshp2022/TGRSHP2022_TechDoc_E.pdf, private roads are classified as mtfcc_code S1740. So we will exclude these roads from the analysis. S1400 & S1500 could also include private roads which may be closed to travel.
wyroads1.head(3)
| permanent_ | source_fea | source_dat | source_d_1 | source_ori | loaddate | interstate | us_route | state_rout | county_rou | federal_la | stco_fipsc | tnmfrc | name | mtfcc_code | intersta_1 | intersta_2 | intersta_3 | us_route_a | us_route_b | us_route_c | state_ro_1 | state_ro_2 | state_ro_3 | SHAPE_Leng | ObjectID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | a631e3a3-278a-49d2-a2c9-f32b8721b8ce | 1106073895730 | {84D2183D-8BB4-4DBB-8848-8FBC90DB5400} | 2017 July MAFTIGER | U.S. Department of Commerce, U.S. Census Burea... | 2017-09-28 | None | None | None | None | None | 56031 | 4 | S Guernsey Rd | S1400 | None | None | None | None | None | None | None | None | None | 0.000469 | 1 | LINESTRING (-104.74206 42.25903, -104.74207 42... |
| 1 | ae22f312-e2b7-4b90-a18b-4f8611272d9a | 1105058887284 | {AD68E48A-6BC8-4301-A05A-3336BCEE0E46} | USFS Roads 04/2021 | U.S. Forest Service | 2021-07-20 | None | None | None | None | None | 56001 | 4 | None | S1400 | None | None | None | None | None | None | None | None | None | 0.004127 | 2 | LINESTRING (-105.47860 42.18840, -105.47854 42... |
| 2 | {5BF6CFD8-FC58-4F2E-9974-2804BCFD2582} | 1105616797687 | {AD68E48A-6BC8-4301-A05A-3336BCEE0E46} | USFS Roads 04/2021 | U.S. Forest Service | 2021-07-20 | None | None | None | None | None | 56005 | 4 | None | S1400 | None | None | None | None | None | None | None | None | None | 0.001327 | 3 | LINESTRING (-105.35928 43.60689, -105.35864 43... |
wyroads1['length'] = wyroads1.geometry.to_crs({'init': 'epsg:5070'}).map(lambda x: x.length)
wyroads1.groupby('mtfcc_code').length.sum().sort_values(ascending=False)
/Users/mhabiger/anaconda3/envs/basic_huntscore/lib/python3.8/site-packages/pyproj/crs/crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string)
mtfcc_code S1400 1.685824e+08 S1500 3.052548e+07 S1200 9.761618e+06 S1740 8.434520e+06 S1100 2.938931e+06 S1630 3.536253e+05 S1640 1.334003e+05 Name: length, dtype: float64
From the above calculation we can see that roads classified as "private" account for 8.4 million meters (~5240 miles) of road. Other possibly private roads are S1400 and S1500 classes. These roads account for the majority of roadway length in Wyoming! For now we will exclude S1500 and S1740 roads in this analysis.
However, from personal experience, I can say that a number of S1400 classified road segments are also private or at least appear private because a sign is posted. Running this notebook with varying mtfcc exclusions will find that Wyoming has potentially between 500,000 to 8 million acres of landlocked land. Determining which segments are public, private or private with a public easement is one of the things the MAPLand Act aims to solve (https://www.congress.gov/bill/117th-congress/house-bill/3113).
wyroads1 = wyroads1[~wyroads1.mtfcc_code.isin(['S1740','S1500'])][['permanent_','geometry']]
wyroads1.columns = ['id','geometry']
wyroads2.head(3)
| permanenti | name | namealtern | trailnumbe | trailnum_1 | sourcefeat | sourcedata | sourceda_1 | sourceorig | loaddate | trailtype | hikerpedes | bicycle | packsaddle | atv | motorcycle | ohvover50i | snowshoe | crosscount | dogsled | snowmobile | nonmotoriz | motorizedw | primarytra | nationaltr | lengthmile | networklen | SHAPE_Leng | ObjectID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 31a64787-fae9-48b3-9edc-2de76087d48b | Mountain Creek Trail | None | None | None | Yellowstone National Park | {A2B9A873-ACAC-4D1C-B014-5B090D31938D} | NPS Trails Update 12/2021 | National Park Service | 2022-01-03 | Terra Trail | Y | N | N | N | N | N | N | N | N | N | None | None | NPS | None | 1.873593 | 3172.813923 | 0.031692 | 1 | LINESTRING (-110.00437 44.32412, -110.00360 44... |
| 1 | 7c1ded41-f5d8-4dc7-8f58-d9fc6c41c64b | Mountain Creek Trail | None | None | None | Yellowstone National Park | {A2B9A873-ACAC-4D1C-B014-5B090D31938D} | NPS Trails 10/2019 | National Park Service | 2019-12-09 | Terra Trail | Y | N | N | N | N | N | N | N | N | N | None | None | NPS | None | 1.436181 | 3172.813923 | 0.022422 | 2 | LINESTRING (-110.00298 44.30561, -110.00326 44... |
| 2 | 57498dbd-a463-49a9-a21f-a0bf29137c8a | Mountain Creek Trail | None | None | None | Yellowstone National Park | {A2B9A873-ACAC-4D1C-B014-5B090D31938D} | NPS Trails 10/2019 | National Park Service | 2019-12-09 | Terra Trail | Y | N | N | N | N | N | N | N | N | N | None | None | NPS | None | 1.315874 | 3172.813923 | 0.025116 | 3 | LINESTRING (-110.08309 44.23790, -110.08282 44... |
wyroads2 = wyroads2[['permanenti','geometry']]
wyroads2.columns = ['id','geometry']
roads = pd.concat([wyroads1,wyroads2])
roads.crs
<Geographic 2D CRS: EPSG:4269> Name: NAD83 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands. - bounds: (167.65, 14.92, -40.73, 86.45) Datum: North American Datum 1983 - Ellipsoid: GRS 1980 - Prime Meridian: Greenwich
We'll need to map some geometries but also want to be able to compare in 2D so we project to AEA
roads['geometry1'] = roads.geometry.to_crs({'init': 'epsg:5070'})
/Users/mhabiger/anaconda3/envs/basic_huntscore/lib/python3.8/site-packages/pyproj/crs/crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string)
First we need to get down to individuals polygons. We explode multipolygon objects
wyland.shape
(19761, 3)
wyland = wyland.explode('geometry')
wyland.shape
(33050, 3)
wyland = wyland.reset_index()
Next, we will generate geohashes to speed up computations. If you aren't familiar with a geohash just google it. It's one of many ways to create a geospatial index that allows partitioning for fast lookup.
wyland['geometry1'] = wyland.geometry
wyland['geometry'] = wyland.geometry.to_crs({'init': 'epsg:4269'})
/Users/mhabiger/anaconda3/envs/basic_huntscore/lib/python3.8/site-packages/pyproj/crs/crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 in_crs_string = _prepare_from_proj_string(in_crs_string)
wyland['geohashes'] = wyland.geometry.map(lambda x: list(polygon_to_geohashes(x, 4, False)))
wyland = wyland.explode('geohashes')
Then we explode on geohashes. The whole purpose is to allow us to compare small chunks of geospatial data rather than the entire set. We don't need to compare a polygon on the Western side of the state to a polygon on the Eastern side of the state and this allows us to compare polygons that are only close.
wyland['acres'] = wyland.geometry1.map(lambda x: x.area*0.000247105)
wyland.shape
(36277, 7)
This index will come in handy later as it allows us to identify a unique polygon. Given there may currently be duplicates but we can handle that in calculations as needed.
wyland['id'] = wyland.level_0.astype(str)+'_'+wyland.level_1.astype(str)
Let's map the 50 largest polygons in the dataset.
osm = folium.Map([42.994,-105.590],zoom_start=6)
mapped = []
n = 0
for i in wyland.sort_values('acres',ascending=False).itertuples():
if i.id not in mapped:
if n<50:
folium.GeoJson(gm.mapping(i.geometry)).add_to(osm)
mapped.append(i.id)
n+=1
osm